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ABSTRACT 

Spherical models of collisionless but quasi-relaxed stellar systems have long 
been studied as a natural framework for the description of globular clusters. Here 
we consider the construction of self-consistent models under the same physical 
conditions, but including explicitly the ingredients that lead to departures from 
spherical symmetry. In particular, we focus on the effects of the tidal field as- 
sociated with the hosting galaxy. We then take a stellar system on a circular 
orbit inside a galaxy represented as a "frozen" external field. The equilibrium 
distribution function is obtained from the one describing the spherical case by 
replacing the energy integral with the relevant Jacobi integral in the presence 
of the external tidal field. Then the construction of the model requires the in- 
vestigation of a singular perturbation problem for an elliptic partial differential 
equation with a free boundary, for which we provide a method of solution to 
any desired order, with explicit solutions to two orders. We outline the relevant 
parameter space, thus opening the way to a systematic study of the properties 
of a two-parameter family of physically justified non-spherical models of quasi- 
relaxed stellar systems. The general method developed here can also be used to 
construct models for which the non-spherical shape is due to internal rotation. 
Eventually, the models will be a useful tool to investigate whether the shapes of 
globular clusters are primarily determined by internal rotation, by external tides, 
or by pressure anisotropy. 

Subject headings: globular clusters: general — methods: analytical — stellar 
dynamics 



Introduction 



Large stellar systems can be studied as collisionless systems, by means of a one-star 
distribution function obeying the combined set of the collisionless Boltzmann equation and 
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the Poisson equation, under the action of the self-consistent mean potential. For elliptical 
galaxies the relevant two-star relaxation times do actually exceed their age; an imprint 
of partial relaxation may be left at the time o f their formation ( if we refer to a p icture 
of formation via incomplete violent relaxation; iLynden-Belll 119671 . Ivan Albadal Il982l ). but 
otherwise they should be thought of as truly collisionless systems, generally characterized 
by an anisotropic pressure tensor. In turn, for globular clusters the relevant relaxation 
times are typically shorter than their age, so that we may argue that for many of them 
the two-star relaxation processes have had enough time to act and to bring them close to 
a thermodynamically relaxed state, with their distribution function close to a Maxwellian. 
This line of argu ments has led to the development of well-known dynamical models for 
globular clusters (lKing|l965l ; Il966l ). 



King models are based on a quasi-Maxwellian isotropic distribution function }k{E) in 
which a truncation prescription, continuous in phase space, is set heuristically to incorporate 
the presence of external tides; but otherwise, they are fully self-consistent (i.e., no external 
fields are actually considered) and perfectly spherical. Empirically, the simplification of 
spherical symmetry is encouraged by the fact that in general globular clusters have round 
appearance. Indeed, as a zero-th order description, these models have had remarkable success 
in applications to observed globular clus ters (e.g., see ISpitzerl 1 198 71 ; iDjorgovski &: Meylan 



1994 ; McLaughlin &; van der Marell 120051 ; and references therein). In recent years, great 
progress has been made in the acquisition of detailed quantitative information about the 
structure of these stellar systems, especially in relation to the measurement of the proper 



motio ns of thousands of individual stars (see Ivan Leeuwen et al.l |2000| ; McLaughlin et al. 



20061 ). with the possibility of getting a direct 5-dimensional view of their phase space. Such 
progress calls for renewed efforts on the side of modeling. More general models would allow 
us to address the issue of the origin of the observed departures from spherical symmetry. 
In fact, it remains to be established which physical ingredient among rotation, pressure 
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As in the case of the study of elliptical galaxies (e.g., see 
references therein), different approaches can be taken to the construction of models. Broadly 
speaking, two complementary paths can be followed. In the first, "descriptive" approach, 
under suitable geometrical (on the intrinsic shape) and dynamical (e.g., on the absence or 
presence of dark matter) hypotheses, the available data for an individual stellar system are 
imposed as constraints to derive the internal orbital structure (distribution function) most 
likely to correspond to the observations. This approach is often carried out in terms of codes 
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that generalize a metho d introduced by ISchwarzschildl (119791 ); for an application to the glob 



ular cluster u Cen, see Ivan de Yen et al.l (120061 ). In the second, "predictive" approach, one 



proposes a formation/evolution scenario in order to identify a physically justified distribution 
function for a wide class of objects, and then proceeds to investigate, by comparison with 
observations of several individual objects, whether the data support the general physical 
picture that has been proposed. Indeed, King models belong to this latter approach. 

The purpose of this paper is to extend the description of quasi-relaxed stellar systems, 
so far basically limited to the spherical King models, to the non-spherical case. There 
are at least three different ways of extending spherical isotropic models of quasi-relaxed 
stellar systems (such as King models), by modifying the distribution function so as to in- 
clude: (i) the explicit presence of a non-spherical tidal field; (ii) the presence of internal 
rotation; (iii) the presence of some pressure anisotropy. As noted earlier in this Introduc- 
tion, these correspond to the physical ingredients that, separately, may be thought to be 
at the origin of the observed non-spherical shapes. We will thus focus on the construction 
of physically justified models, as an extension of the King models, in the presence of ex- 
ternal tides and, briefly, on the extension of the models to the presence of rigid internal 
rotation. A firs t-order analysis of the tr iaxial tidal problem addressed in this paper was car- 



ried out by lHeggie Ramamanil (11995) and the e ffect of a "frozen" tidal field on (initially) 



spherical King models was studied by IWeinberd (119931 ) using N-body simulations. For the 
axisymmetric problem associated with the presence of internal rotation, a first-order analysis 
of a simply truncate d Max wellian distribution perturbed by a rigid rotation was given by 
Kormendv fc Anandl (119711); different mo dels whe re differential rotation is considered were 
proposed by iPrendergast fc Tomerl (119701 ) and by IWilsonl (119751 ). also in view of extensions 
to the presence of pressure anisotropy, which goes beyond the scope of this paper. Models 
that represent a direct generalization of the King family to the case with differential rota- 
tion have also been examine d , with particular attention t o their thermodynamic properties 
(ILagoute fc Longarettilll996l ; lLongaretti fc Lagoutdll996l ). In principle, the method of so- 
lution that we will present below can deal with the extension o f other spherical isotropic 
models with fi nite size that are not of the King fo rm (e.g., see Iwoolley fc Dickens 1962 
Davoust 1977 ; see also the interesting sugg estion by lMadsenlll996l l 



The study of self-consistent collisionless equilibrium mo dels ha s a long tradition not only 
in st ellar dynamics, but also in plasma physics (e.g., see iHarrid Il962l ; lAttico fc Pegoraro 
19991 ). We note that in both research study in the presence of external fields, 

especially when the external field is bound to break the natural symmetry associated with 
the one-component problem, is only rarely considered. 



The paper is organized as follows. Section 2 introduces the reference physical model, in 
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which a globular cluster is imagined to move on a circular orbit inside a host galaxy treated 
as a frozen background field; the modified distribution function for such a cluster is then 
identified and the relevant parameter space defined. In Sect. 3 we set the mathematical 
problem associated with the construction of the related self-consistent models. For models 
generated by the spherical fx(E), Sect. 4 and Appendix A give the complete solution in terms 
of matched asymptotic expansions. Alternative methods of solutions are briefly discussed in 
Sect. 5. The concluding Section 6 gives a summary of the paper, with a short discussion of 
the results obtained. In Appendix B we show how the method developed in this paper can 
be applied to construct quasi-relaxed models flattened by rotation in the absence of external 
tides. In Appendix C we show how the method can be applied to other isotropic truncated 
models, different from King models. 

Technically, the mathematical problem of a singular perturbation with a free boundary 
that is faced here is ver y similar to the problem noted in the theory of rotating stars , 



starting with Milne (see iTassoull Il978l: iMilnel Il923l: Ichandrasekhar! 1 19331 : iKrogdahll If942 



Chandrasekhar fc Lebovita 19621 ; Monaghan fc Roxburgh 1965 ). The problem was initially 



dealt with inadequate t ools; a satisfactory solution of the singular perturbation problem was 
obtained only later, by ISmithl (119751 ; Il976l ). 



2. The physical model 

2.1. The tidal potential 

As a reference case, we consider an idealized model in which the center of mass of a 
globular cluster is imagined to move on a circular orbit of radius Ro, characterized by orbital 
frequency Q, inside a host galaxy. For simplicity, we focus on the motion of the stars inside 
the globular cluster and model the galaxy, taken to have very large mass, by means of a 
frozen gravitational field (which we will call the galactic field, described by the potential 
with a given overall symmetry. This choice makes us ignore interesting effects that 
are generally present in the full interaction between a "satellite" and a galaxy; in a sense, 
we are taking a complementary view of an extremely complex dynamical situation, with 
respect to other investigations, such as those that lead to a discussion of the mechanism 
of dynamical friction (in which the globular cluster or satellite is modeled as a rigid body 



and the stars of the galaxy are taken as the " live" component; see iBontekoe fc van Albada 



19871 ; iBertin et al.ll2003l ; lArena fc Bertinll2007l ; and references therein). Therefore, we will be 



initially following the picture of a restricted three-body problem, with one important difference, 
that the "secondary" is not treated as a point mass but as a "live" stellar system, described 
by the cluster mean- field potential $c- hi this extremely simple orbital choice for the cluster 
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center of mass, in the corotating frame the associated tidal field is time-independent and so 
we can proceed to the construction of a stationary dynamical model. 

We consider the galactic potential $g to be spherically symmetric, i.e. &g — ®g(R), 
with R = \JX 2 + Y 2 + Z 2 , in terms of a standard set of Cartesian coordinates (X, Y, Z) , so 
that Q 2 = (d§G(R)/dR)n / Ro- Let (X,Y) be the orbit plane of the center of mass of the 
cluster. We then introduce a local rotating frame of reference, so that the position vector 
is given by r = (x,y,z), with origin in the center of mass of the cluster and for which the 
x-axis points away from the center of the galactic field, the y-axis follows the direction of 
the cluster rotation in its orbit around the galaxy, and the z-axis is perpendicular to the 
orbit plane (according to the right-hand rule). In such rotating local fra me, the relevant 



Lagr angian, describing the motion of a star belonging to the cluster, is (cf. IChandrasekhar 



1941 Eq. [5.510]): 

C = ^{x 2 + y 2 + z 2 + n 2 [(R + x) 2 + y 2 ] + 2{l(R + x)y - 2Slxy} - $ G (R) - $ c (x, y, z) 



where R = a/ (-Ro + x) 2 + y 2 + z 2 and the terms responsible for centrifugal and Coriolis 
forces are explicitly displayed. 

If we suppose that the size of the cluster is small compared to Ro, we can adequately 
represent the galactic field by its linear approximation with respect to the local coordinates 
(the so-called "tidal approximation"). The corresponding equations of the motion for a single 
star in the rotating local frame are: 

x - 2Vty - (An 2 - k 2 )x = , (2) 

y + 2Qx = -—^, (3) 
dy 

z + n 2 z = ^, (4) 

oz 

where k is the epicyclic frequency at Rq, given by n 2 = 3Q 2 + (gP$ g / 'dR 2 )^. Note that the 
assumed symmetry for $c introduces a cancellation between the kinematic term yQ 2 and 
the gradient of the galactic potential d^ajdy and makes the vertical acceleration —d$>a/dz 
approximately equal to — zfl 2 . 

These equations admit an energy (isolating) integral of the motion, known as the Jacobi 
integral: 

H = ^(x 2 + y 2 + z 2 ) + $ T + $ c , (5) 

where 

$ T = -ti 2 (z 2 - vx 2 ) (6) 
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is the tidal potential. Here v = 4 — k 2 /Q 2 is a generally positive dimensionless coefficient. 



Thus, at the level of single star orbits, we note that, in general, the tidal potential leads 
to a compression in the z-direction, a stretching in the x-direction, and leaves the y-direction 
untouched. The tidal potential is static, breaks the spherical symmetry, but is characterized 
by reflection symmetry with respect to the three natural coordinate planes; strictly speaking, 
the symmetry w ith respect to (y, z) is applicable only in the limit of an infinitely massive 
host galaxy (see ISpitzerlll987l ). In turn, we will see that the geometry of the tidal potential 
induces a non-spherical distortion of the cluster shape collectively, in particular an elongation 
along the x-axis and a compression along the z-axis. In practice, the numerical coefficient 
v that defines quantitatively the induced distortion depends on the galactic potential. We 
recall that we have v = 3 for a Keplerian potential, v = 2 for a logarithmic potential, while 
for a Plummer model the dimensionless coefficient depends on the location of the circular 
orbit with respect to the model scale radius b, with is(Rq) = 3Rl/(b 2 + Rq) (for a definition 
of the Plummer model see, e.g.. iBertinl 120001 ). 



Different assumptions on the geometry of the galactic field can be treated with tools 
similar to those developed here, leading to a similar structure of the equations of the motion, 
with a slight modification of the tidal field. In particular, for an axisy mmetric galactic field , 



the tidal potential d iffers from the one obtained here only by the z-term (jChandrasekharlll942 
Heggie &: Hutl 120031 ). This case is often considered, for example by r eferring to a globular 



clust e r in circular orbit on the (axisymmetric) disk of our Galaxy (see lHeggie fc Ramamani 
19951 ; lErnst et al.ll2008l ). for which $^ is then formulated in terms of the Oort constants. 



In the physical model outlined in this Section, the typical dynamical time associated with 
the star orbits inside the cluster is much smaller than the (external) orbital time associated 
with Q. Therefore, in an asymptotic sense, the equilibrium configurations that we will 
construct in the rest of the paper can actually be generalized, with due qualifications, to 
more general orbits of the cluster inside a galaxy, provided we interpret the results that we 
are going to obtain as applicable only to a small piece of the cluster orbit. 



2.2. The distribution function 

As outlined in the Introduction, we wish to extend the description of quasi-relaxed 
stellar systems (so far basically limited to spherical models associated with distribution 
functions / = f(E), dependent only on the single-star specific energy E = t> 2 /2 + <3>c) to the 
non spherical case, by including the presence of a non-spherical tidal field explicitly. Given 
the success of the spherical King models in the study of globular clusters, we will focus 
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on the extension of models based on fx(E), which is defined as a "lowered" Maxwellian, 
continuous in phase space, with an energy cut-off which implies the existence of a boundary 
at the truncation radius r tr . 

Therefore, we will consider (partially) self-consistent models characterized by the dis- 
tribution function: 

/* (if) = A[exp(-aH) - exp(-aff )] (7) 

if if < if and fx(H) = otherwise, in terms of the Jacobi integral defined by Eq. ((SJ). Here 
ffo is the cut-off value for the Jacobi integral, while A and a are positive constants. 

In velocity space, the inequality H < H identifies a spherical region given by < v 2 < 
2ip(r)/a, where: 

V(r) = a{H - [$ c (r) + $ T (x, z)}} (8) 

is the dimensionless escape energy. Therefore, the boundary of the cluster is defined as the 
relevant zero- velocity surface by the condition ^(r) = and is given only implicitly by an 
equipotential (Hill) surface for the total potential $c + ^r; in fact, its geometry depends on 
the properties of the tidal potential (of known characteristics; see Eq. [6]) and of the cluster 
potential (unknown a priori, to be determined as the solution of the associated Poisson 
equation). 

The value of the cut-off potential ifo should be chosen in such a way that the surface 
that defines the boundary is closed. The last (i.e., outermost) closed Hill surface is a critical 
surface, because it contains two saddle points that represent the Lagrangian points of the 
restricted three-body problem outlined in the previous subsection. From Eqs. (J2D- (SJ) , we 
see that such Lagrangian points are located symmetrically with respect to the origin of the 
local frame of reference and lie on the x-axis. Their distance from the origin is called the 
tidal radius, which we denote by r>r, and can be determined from the condition: 

^(r T ,0,0)=0. (9) 

If, as a zero-th order approximation, we use a simple Keplerian potential for the cluster 
potential $c, we recover the classical expression (e.g., see Spitzer 1987): 

• (10) 

where M is the total mass of the cluster. 

As for the spherical King model, the density profile associated with the distribution 
function ([7]) is given by: 

p(i/f) = Ae*j =Ap(4>) , (11) 
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where A = 87iA\/2e~ aHo / (3a 3 ^ 2 ). We recall that the incomplete gamma function has non- 
negative real value only in correspondence to a non- negative argument. In the following, we 
will denote the central density of the cluster by po — Apfi?), where \1/ = ^(O) is the depth 
of the central potential well. 



2.3. The parameter space 

The models defined by fx{H) are characterized by two physical scales (e.g., the two 
free constants A and a, or, correspondingly, the total mass M and the central density po of 
the cluster) and two dimensionless parameters. The first dimensionless parameter can be 
defined, as in the spherical King models, to measure the concentration of the cluster. We 
can thus consider the quantity introduced at the end of the previous subsection, or we 
may refer to the commonly used concentration parameter: 

C = log(7v/r ) , (12) 

where r = A/9/(47rGpo a ) is a scale length related to the size of the core and r tr is the 
truncation radius of the spherical King model associated with the sam e value of th e central 



potential well \I> (the relation between C and \l/ is one-to-one; e.g., see iBertiru l2000l ) 



The second dimensionless parameter characterizes the strength of the (external) tidal 

field: 

The definition arises naturally from the dimensionless formulation of the Poisson equation 
that describes the (partially) self-consistent problem (to be addressed in the next Section). 

In principle, for a given choice of the dimensional scales (A and a) the truncation radius 
or the concentration parameter of a spherical King model can be set arbitrarily. In practice, 
the physical motivation of the models suggests that the truncation radius r tr should be 
taken to be of the order of (and not exceed) the tidal radius Pr, introduced in the previous 
subsection (see Eq. [5]). We may thus define an extension parameter, as the ratio between 
the truncation radius of the corresponding spherical model and the tidal radius r^: 

5=^. (14) 
r T 

For a given value of the central potential well there exists a (maximum) critical value for 
the tidal strength parameter, which we will denote by e cr , corresponding to the maximum 
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value for the extension parameter 5 cr , which can be found by solving the system: 

^,0,0;^) =0 (15) 
7/>(r T ,0,0;e cr 



. 



From this syst em, if we use the zero-th order Keplerian approximation for <3>(7, we find that 



2/3 (see ISpitzerl Il987h . 



For our two-parameter family of models we thus expect two tidal regimes to exist. For 
models characterized by the pairs e) near the critical condition 5 ~ 5 cr the tidal distortion 
should be maximal, while for models with pairs well below criticality only small departures 
from spherical symmetry should occur. A thorough exploration of the parameter space will 
be carried out in a separate paper ( Varri fc Bertin, in preparatio n). In closing, we note that 
the models proposed and studied by iHeggie fc Ramamanil (119951 ) correspond to the pairs in 
parameter space that we have called critical. 



3. The mathematical problem 

The (partially) self-consistent models associated with the distribution function defined 
by Eq. (jTJ) are constructed by solving the relevant Poisson equation. In terms of the dimen- 
sionless escape energy ip, given by Eq. (jSj), the Poisson equation (for ip > 0) can be written 
as: 

72 /., i „>k \ 9 p 9 p(4>) 



r o Po 



(16) 



where r is the scale length introduced in Sect. 12.31 (see Eq. [12]). We then rescale the 
coordinates and introduce the dimensionless position vector r = r/ro, so that V 2 = TqV 2 
and a$r = eT = 9e(z 2 — ux 2 )/2, where we have made use of the tidal parameter introduced 
in Eq. ([TBI . Therefore, the Poisson equation, for ip > 0, can be written in dimensionless 
form as: 



-9 



while for negative values of ip we should refer to: 

V 2 ?/> = -9e(l 
i.e. the Laplace equation V 2 (a$c') = 0. 



+ c(l 



(17) 



(18) 



The mathematical problem is completed by specifying the appropriate boundary condi- 
tions. As for the spherical King models, we require regularity of the solution at the origin 



^(0) = * 



(19) 
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and, at large radii: 



which corresponds to a$c* ~~ * 0. 



W(0) = , (20) 
^ + eT -> a# , (21) 



Poisson and Laplace domains are thus separated by the surface defined by tp = which 
is unknown a priori] in other words, we have to solve an elliptic partial differential equation 
in a free boundary problem. 

In the ordinary differential problem that characterizes the construction of spherical 
models with finite mass, the condition of vanishing cluster potential at large radii (together 
with the regularity conditions at the origin) overde termines the problem, w hich can then be 



seen as an eigenvalue problem (e.g., see Sect. 2.5 in lBertin fc Stiaveih1ll993l ). Indeed, for the 
King models the integration of the Poisson equation from the origin outwards, with "initial 
conditions" ffT^ - fpZUl) . sets the relation between the ratio r tr /ro and \1/ in order to meet the 
requirement (l2Tj) . with r tr /ro thus playing the role of an "eigenvalue". 

In the more complex, three-dimensional situation that we are facing here, the existence 
of two different domains, internal (Poisson) and external (Laplace), suggests the use of the 
method of matched asymptotic expansions in order to obtain a uniform solution across the 
separation free surface. The solution in the internal and external domains are expressed as 
an asymptotic series with respect to the tidal parameter e, which is assumed to be small 
(following the physical model described in the previous Section): 

^(r;e)=J2y4 nt \r)e k , (22) 



k=0 



oo _ 

rt^E^) 6 *. ( 23 ) 

fc=0 

with spherical symmetry assumed for the zero-th order terms. The internal solution should 
obey the boundary conditions (|T9j) - (j20"I) . while the external solution should satisfy Eq. ( l2~TT) . 
The two representations should be properly connected at the surface of the cluster. 

On the other hand, for any small but finite value of e the boundary, defined by ip = 0, will 
be different from the unperturbed boundary, defined by ipo — 0, so that, for each of the two 
representations given above, there will be a small region in the vicinity of the surface of the 
cluster where the leading term is vanishingly small and actually smaller than the remaining 
terms of the formal series. Therefore, we expect the validity of the expansion to break down 
where the second term becomes comparable to the first, i.e where ipo = O(e). This region 
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can be considered as a boundary layer, which should be examined in "microscopic" detail 
by a suitable rescaling of the spatial coordinates and for which an adequate solution tp( lay \ 
expressed as a different asymptotic series, should be constructed. To obtain a uniformly valid 
solution over the entire space, an asymptotic matching is performed between the pairs (ip^ nt \ 
^(lay)^ anc j ^(iay) ^ ^(ext)^ th us leading to a solution ip, obeying all the desired boundary 
conditions, in terms of three different, but m atched, representations. This method of solution 
is basically the same method proposed by ISmithl (119751 ) for the analogous mathematical 
problem that arises in the determination of the structure of rigidly rotating fluid polytropes. 



4. Solution in terms of matched asymptotic expansions 



The complete solution to two significant orders in the tidal parameter is now presented. 
The formal solution to three orders is also displayed becau se of the requirements of the Van 
Dyke principle of asymptotic matching (cf. IVan Dykdll975l . Eq. [5.24]) that we have adopted. 



4.1. Internal region 

If we insert the series ( )22|) in the Poisson equation (TlTl) . under the conditions (|T9|) - (f20|) . 
we obtain an (infinite) set of Cauchy problems for ip^. The problem for the zero-th order 
term (i.e., the unperturbed problem with e = 0) is the one defining the construction of the 
spherical and fully self-consistent King models: 

p (ib {int) ) 

4 nt) " + l4 intY = -V-^^J 1 > ( 24 ) 

with ipQ nt \o) = ^ and ip^ 71 ^ (0) = 0, where the symbol ' denotes derivative with respect to 
the argument f. We recall that the truncation radius f tr , which defines the boundary of the 
spherical models, is given implicitly by ipQ nt \f tr ) = 0. 



Let us introduce the quantities: 

9 d?p 



p(^) dipi 



(25) 



These quantities depend on f implicitly, through the function x/jq \ in turn, the dependence 
on ^ is both explicit (through the term p(\l/)) and implicit (because the function ipQ nt \r) 
depends on the value of \l/). For convenience, we give the expression of the first terms of the 
sequence (cf. Eq. [11]): R x = [9/p(*)] + (^fV 2 ], ^2 = ^ + 27^)^ /[2p^)l 
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R 3 = R 2 + 27(4 in ' ) )" 1/2 /[4p(^)]. Note that for f[ !nl) -> 0, i.e. for f -> f tr , R x -> 0, i? 2 -> 0, 
while for j > 3 the quantities .R.,- actually diverge. This is one more indication of the singular 
character of our perturbation analysis, which brings in some fractional power dependence on 
the perturbation parameter e (see also expansion [46] ). 

Therefore, the equations governing the next two orders (for with k — 1, 2) can be 
written as: 

_ V 2 + i?!(r; ^ int) = -9(1 - v) (26) 

V 2 + R x {r- *)] * 3 = -i? 2 (r; *)(^f nt) ) 2 (27) 

with ^i m4) (0) = ^? n * } (0) = and V^J int) (0) = V^ n ' } (0) = 0. The equation for k = 3 is 
recorded in Appendix A.l, where we also describe the structure of the general equation for 

, (int) 

n ■ 

For any given order of the expansion, the operator acting on the function ip^ (see the 
left-hand side of Eqs. [2B] and [27]) is the same, i.e. a Laplacian "shifted" by the function 
Ri(f; \&). If we thus expand every term ipk(r) in spherical harmonics^: 



1=0 m=-l 



the three-dimensional differential problem is reduced to a one-dimensional (radial) problem, 
characterized by the following second order, linear ordinary differential operator: 



d 2 2 d 1(1 + 1) T 
V l = — + ^-- + R x {r- * . (29) 



In general, for a fixed value of I, two independent solutions to the homogeneous problem 
T>if = are expected § to behave like r l and l/f l+1 for f — > 0. Because of the presence of 
i?i(f; \1/), solutions to equations where X>z appears have to be obtained numerically. 

For k = 1 (see Eq. [26]) we thus have to address the following problem. For I = 0, the 
relevant equation is: 

Z>o/oo = -9(1-1/), (30) 



■""We use orthonormalized real spherical harmonics with Condon-Shortley phase; with respect to the 
toroidal angle </>, they are even for m > and odd otherwise. 

2 We note that i?i(0, = 9[1 + * 3 / 2 /pWL i.e. a numerical positive constant. Therefore, for r — > 
the operator 2?/ tend s to the operator associated with the spherical Bessel functions of the first and of the 
second kind (e.g., see Abramowitz fc Stegunl 19641 Eq. [10.1.1] for the equation and Eqs. [10.1.4] and [10.1.5] 
for the limiting values of the functions for small argument). 
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where /oo = V100 / V^r , with /oo(0) = /o (0) = 0. Here we do not have to worry about 
including solutions to the associated homogeneous problem, because one of the two indepen- 
dent solutions would be singular at the origin and the other would be forced to vanish by 
the required condition at f = 0. For I > 1 we have: 

Wffi = (31) 

with ^i*]^(0) = T^^Im (0) = 0- Both Eq. ( J3il and the associated boundary conditions are 
homogeneous. Therefore, the solution is undetermined by an m-dependent multiplicative 
constant: V'umO^) = A i?n 7;(r), with 7/(f) ~ f l for f — > (the singular solution is excluded 
by the boundary conditions at the origin). Then the complete formal solution is: 

oo I 

^ nt) H) = /oo(r) + E E ^7*(r)yU0, 0) , (32) 

1=1 m=—l 

where the constants are ready to be determined by means of the asymptotic matching with 
i/jf ay \r) at the boundary layer. 



For k = 2 (see Eq. [27] ) the relevant equations are: 



(int) 2 



I in 

(int) . 



(33) 



where on the right-hand side the function ip[ is that derived from the solution of the first 
order problem (which shows the progressive character of this method for the construction 
of solutions). In Appendix A. 2 the equations for the six relevant harmonics are displayed 
explicitly. The boundary conditions to be imposed at the origin are again homogeneous: 
V^lm vO) = V4*Im (0) = 0- F° r a fixed harmonic (l,m) with I > 0, the general solution 
of Eq. fl33|) is the sum of a particular solution (which we will denote by gi m (r)) and of a 
regular solution to the associated homogeneous problem given by Eq. (|3T|) (which we will call 
Bim)i{r), with 7;(f) the same functions introduced for the first order problem). Obviously, 
the particular solution exists only when Eq. (I3"3l is non-homogeneous, i.e. only for those 
values of (I, m) that correspond to a non- vanishing coefficient in the expansion of (4 mt) ) 2 m 
spherical harmonics. As noted in the first order problem (k = 1), the associated homogeneous 
problem for / = has no non-trivial solution. Then we can express the complete solution as: 



4 nt \i) = 9oo(r) + E E ^im{r) + B lmll (r)]Y lm (9, 0) , (34) 

i=l m=—l 



where B[ m are constants to be determined from the matching with the boundary layer. 
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Similarly, for k — 3 the solution can be written as: 

oo I 

4 nt \r) = h 00 (r) + £ £ ^(f) + C imli(rWi m (0, <P) , (35) 

1=1 m=—l 

where hi m are particular solutions and Ci m are constants, again to be determined from the 
matching with the boundary layer. 

Because the differential operator T>i and the boundary conditions at the origin are the 
same for the reduced radial problem of every order, we have thus obtained the general 
structure of the solution for the internal region (see Appendix A.l). 



4.2. External region 

Here we first present the general solution and then proceed to set up the asymptotic 
series f[2"3"j) . 

The solution to Eq. (|T8l) describing the external region, i.e. in the Laplace domain, can 
be expressed as the sum of a particular solution (— eT(f)) and of the solutions to the radial 
part of the Laplacian operator consistent with the boundary condition (I2TI) : 

\ °° 1 R 

^(?) = « - £ - £ E ^Ue, <P) - eT{v) . (36) 

1=1 m=—l 

Here we note that the tidal potential contributes only with spherical harmonics of order 
/ = 0, 2 with even values of m: 

T 00 (f) = -3^(u-l)f\ (37) 
T 20 (r) = 3J|(2 + v)f 2 , (38) 

T 22 (r) = -3^z/f 2 . (39) 

At this point we can proceed to set up the asymptotic series, by expanding the constant 
coefficients a, A, and $ m with respect to e: 

a = aH Q = a + a x e + ^ya 2 e 2 + •• , (40) 
A = A + A 1 e+^A 2 e 2 + .. , (41) 
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0i 



1, 2 1 3 
Im. — a lm£ + ^jO/ m 6 + ~^ C lm e + 



(42) 



The last expansion starts with a first order term because the density distribution of the 
unperturbed problem is spherically symmetric. 



For convenience, we give the explicit expression of the external solution up to third 



order: 




Too(f) 
2^ 



i 



+ 



2! 



A 2 

"2 ; 



EE 

1=1 m=—l 



Y 



lm\ 



6 2 + 



EE 

1=1 m=—l 
1 



dim 
rl+1 



+ T, m (f 



Y 



lm \ 



" .; 



(=1 m=—l 



4.3. Boundary layer 

The boundary layer is the region where the function ip becomes vanishingly small. Since 
the unperturbed gravitational field at the truncation radius is finite, ip' {f tr ) ^ 0, for any 
value of based on a Taylor expansion of ipo about f = f tr we may argue that the region in 
which the series ( 1221) breaks down can be defined by f tr — f = 0(e). In this boundary layer 
we thus introduce a suitable change of variables: 

77=^, (44) 
e 

take the ordering tp( lay ^ = (9(e), and thus rescale the solution by introducing the function 
t = ip( lay '/e. For positive values of r the Poisson equation ( flTl) thus becomes: 

tt^ - tj + 77 r = -— — ep(er) - 9e 2 l - i/) , 45 

Of] rtr — or] (r tr — erj) 2 p(W) 

where A 2 is the angular part of the Laplacian in spherical coordinates. For negative values of 
t we can write a similar equation, corresponding to Eq. ffl8l) . which is obtained from Eq. (j4"5l) 
by dropping the term proportional to p(er). 

With the hel p of the asymptotic exp ansion for small argument of the incomplete gamma 



function (e.g., see lBenderfe Orszadll999l . Eq. [6.2.5]), we find: 



p(er)~ 2 r 5/ 2e 5/ 2 + | r 7/ 2e 7/ 2 + __ (46) 

so that, within the boundary layer, the contribution of p(er) (which is the one that dis- 
tinguishes the Poisson from the Laplace regime) becomes significant only beyond the tidal 
term, clS db correction 0(e 7 ' 2 ). 
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Therefore, up to 0(e 2 ) we can write: 



T = T + T X e 



2! 



r 2 e 2 . 



(47) 



To this order, which is required for a full solution up to k = 2 of the global problem (see 
Eqs. [22] and [23]), by equating in Eq. (|45|) the first powers of e separately, we obtain the 
relevant equations for the first three terms: 

d 2 T 



drf 

d 2 n 2dro 
dr] 2 f-tr dr] 



d 2 T 2 

dr] 2 



4 

r tr 



9ti r] dr 
dr] r tr dr] 



- ttA 2 t - 18(1 - v) 



tr 



The equations are easily integrated in the variable rj, to obtain the solutions: 

T = F {6,<f>)ri + G (e,<i>) , 



Tl 



2F n 



T-2 



-V 



1 A 2 F (e,<p)v 3 + 2lli 



tr 



Sf 2 



Ttr 



-v 



-9(1 - v)rf 



-A 2 G (6,<j ) )r] 2 + F 2 (9,<P)r] + G 2 



(48) 

(49) 
(50) 

(51) 
(52) 



(53) 



The six free angular functions that appear in the formal solutions will be determined by the 
matching procedure. 



4.4. Asymptotic matching to two orders 

In order to obtain the solution, we must perform separately the relevant matching for 
the pairs (i/)( tnt \ ip( lay )) and (ip( lay \ip( ext ' ) ). We follow the Van Dyke matching principle, 
which requires that we compare the second order expansion of the internal and external 
solutions with the third order expansion of the boundary layer solution. The full procedure 
is described in Appendix A. 3. 

To first order (i.e., up to k = 1 in series [22J and [23]), from the matching of the pair 
{ip^ nt \ ip(! a v>) we find the free angular functions of (1511) and (1521) : 

F o (0,<t>) = -4 int) '(r tr ) , (54) 
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G (6,<p) = ^ nt \rtr,e,<P) , (55) 

Fi{eA) = - J ^-{hr,e,<i>) , (56) 
G x {9A) = \4 nt \hr,e,4>) . (57) 

From the matching of the pair (ip( ext \ ^C a ^)) we connect ip( ext *> to the same angular functions, 
thus proving that the matching to first order is equivalent to imposing continuity of the 
solution up to second order and of the first radial derivative up to first order. This allows us 
to determine the free constants that are present in the first two terms of (1431) and in (|32|) : 

a = ^- , (58) 

r tr 

Ao = rl4 inty (r tr ) , (59) 
«i = /oo(r tr ) + r tr f 00 (r tr ) + , (60) 

2V7T 

Ai = r tr f 00 {r tr ) + = , (61) 



7T 

, 5T 2m (f ir ) 

= -^(^) + 3-, 2 (f ir ) • (62) 

a 2m = -rt r [A 2m j 2 {hr) + T 2m {hr)) ■ (63) 

Note that Ai m = ai m = if / ^ 2, for every value of m, and that the constants for / = 2 are 
non- vanishing only for m = 0, 2. The constants that identify the solution are thus expressed 
in terms of the values of the unperturbed field ip^ 1 ^ , of the "driving" tidal potential TJ m , 
and of the solutions / o an d 72 (see Eqs. [SD] and [HI]) taken at f = f tr . 

The boundary surface of the first order model is defined implicitly by i/j^ xt ^ {r)+ip^ xt " > (f, 8, < 
0, i.e. the spherical shape of the King model is modified by monopole and quadrupole con- 
tributions, which are even with respect to toroidal and poloidal angles and characterized 
by reflection symmetry with respect to the three natural coordinates planes. As might have 
been expected from the physical model, the spherical shape is thus modified only by spherical 
harmonics (/, m) for which the tidal potential has non- vanishing coefficients. Mathematically, 
this is non-trivial, because the first order equation in the internal region Eq. (126|) is non- 
homogeneous only for Z = 0; the quadrupole contribution to the internal solution is formally 
"hidden" by the use of the function if) (which includes the tidal potential) and is unveiled 
by the matching which demonstrates that A 2m with m = 0, 2 are non- vanishing. 

The first order solution can be inserted into the right-hand side of Eq. (133]) to gen- 
erate non-homogeneous equations (and thus particular solutions) only for I = 0,2,4 and 
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corresponding positive and even values of m (see Appendix A. 2). We can thus proceed 
to construct the second order solution in the same way described above for the first order 
solution. From the matching of the pair (ip( int \ ip( lay ty we determine the missing angular 
functions: 

F 2 (e,<p) = -^-(f tr ,e, ( j ) ) , (64) 

G 2 (M) = ^f%r,M) , (65) 

which are then connected to the properties of ijj( ext ^ by the matching of the pair 
This is equivalent to imposing continuity of the solution up to third order and of the first 
radial derivative up to second order and leads to the determination of the free constants that 
appear in the third term of (|43l) and in (134"]) : 

«2 = 9oo(ftr) + rtr9oo(hr) , (66) 

^2 = f 2 tr g' m {f tr ) , (67) 
S2 = r tr g' 2m (f tr ) + 3g 2m {f tr ) ^ 

hm = -rl r [g2m{r t r) + B 2m <y 2 (r tr )} , (69) 
£ 4 __ r tr g' im {r tr ) + 5g 4m (fy) 

&4m = -ft r [g4m(ftr) + B Am ^ A (r tr )] . (71) 

Here Bi m = b[ m = if I ^ 2,4 for every value of m; the only non-vanishing constants with 
I = 2, 4 are those with even m. 

Therefore, the second order solution has non-vanishing contributions only for I = 0,2, 4, 
i.e. for those harmonics for which the particular solution to Eq. (127)) is non-trivial. By 
induction, it can be proved (see Appendix A. 4) that the k-th order solution is characterized 
by I = 0,2,. .,2k harmonics with corrisponding positive and even values of m. In reality, 
the discussion of the matching to higher orders (k > 3) would require a re-definition of the 
boundary layer, because the density contribution on the right-hand side of Eq. ()45)) (for 
positive values of r) comes into play. The asymptotic matching procedure carries through 
also in this more complex case but, for simplicity, is omitted here. We should also keep in 
mind that in an asymptotic analysis the inclusion of higher order terms does not necessarily 
lead to better accuracy in the solution; the optimal truncation in the asymptotic series 
depends on the value of the expansion parameter (in this case, on the value of e) and has to 
be judged empirically. 
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In conclusion, starting from a given value of the King concentration parameter \1/ and 
from a given strength of the tidal field e, the uniform triaxial solution is constructed by nu- 
merically integrating Eqs. ($3$, (13"U1) . (13"T1) . and (13"3"j) and by applying the constants derived in 
this subsection to the asymptotic series expansion (1221) - (1251) . The numerical integrations can 
be performed efficiently by means of standard Runge-Kutta routines. The boundary surface 
of the model is thus defined by ^ xt \f) + ip[ ext \f, 9, 0)e + -?/4 e ^(T> 0) e V2 = 0, while the in- 
ternal density distribution is given by p = p(ipQ int \r)+il)[ mt \r , 9, 0)e+V>2 m ^(r, 9, 0) e V2), with 
the function p defined by Eq. (TTTT) . Any other "observable" quantity can be reconstructed by 
suitable integration in phase space of the distribution function fx^H) defined by Eq. ©, with 
H defined by Eq. ©, and $ T + $ c = H - [^ int \r) + rf nt) (f , 9, (j))e + ^ nt \r, 9, (j))e 2 /2]/a. 

In Fig. 1 and Fig. 2 we illustrate the main characteristics of one triaxial model con- 
structed with the method described in this Section. 



5. Alternative methods of solution 

5.1. The method of strained coordinates 

The mathematical problem described in Sect. 3 can also be solved by the method of 
strained coordinates, an al ternative metho d usually applied to non-linear hyperbolic dif fer- 



ential equations (e.g., see IVan Dykd Il975l . Chapter 6) and considered by ISmith Jl976h 



in 

the solution of the singular free-boundary perturbation problem that arises in the study of 
rotating polytropes. 

Starting from a series representation of the form (|22|) and (|23|) for the solution defined in 
the Poisson and Laplace domains, respectively, a transformation is considered from spherical 
coordinates (f,9,4>) to "strained coordinates" (s,p,q): 

f = s + eri(s,p,g) + ^e 2 f 2 (s,p,q) + ... 
9 = p 

<t> = q, (72) 

where fk(s,p, q) are initially unspecified straining functions. We note that the zero-th order 
problem is defined by the same Eq. (124)) with the same boundary conditions but with the 
variable f replaced by s. The unperturbed spherical boundary in the strained space is defined 
by s = sq, where ^c) ( s o) — 0. To each order, the effective boundary of the perturbed 
configuration remains described by the surface s = sq, while in physical coordinates the 
truncation radius actually changes as a result of the straining functions that are determined 
progressively. 



-20- 



The Laplacian expressed in the new coordinates, V 2 , can be written as an asymptotic 
series: V 2 = L + eLi + l/2e 2 L 2 + where Lf. are linear second order operators B in which 
fj(s,p,q) (with j = l,..,k) and their derivatives appear. For convenience, we record the 
zero-th and first order operators: 

d 2 2 d 

L ^d^ + STs' (73) 

t dri\ d 2 ( d 2 f 1 t 2 dfx i 1 2 „ i 2 A \ d 



^ 1 \ ds ) ds 2 \ ds 2 s ds s 2 ^ ^ s 2 ^ 1 ) ds ' 

where A 2 is the standard angular part of the Laplacian, written with angular coordinates 
{p,q). The general fc-th order operator can be decomposed as L k = L\ + Fk, where F^ is, 
in turn, a second order operator in which fj(s,p, q) (with j = 1, .., k — 1) appear and L\ is 
defined as in Eq. (I73|) but with f k (s,p, q) instead of fi(s,p, q); these operators appear in the 
relevant equation for ip^ nt ^: 

[L Q + R x {^e)]4 nt) = L k ^ (75) 

which corresponds to the general fc-th order equation of the previous method. 

Following a set of constraints that guarantee the regularity of the series (1221) in the 
strained space, the equations that uniquely identify the straining functions to any desired or- 
der can be found and solved numerically; structurally, they somewhat correspond to Eqs. (T3TT) 
and (1551) in Sect. 4.1. Therefore, the internal and external solutions can be worked out and 
patched by requiring continuity of the solution and of the first derivative with respect to the 
variable s at the boundary surface defined by s = sq, in general qualitative analogy with the 
method described in the previous section. 

This method is formally more elegant than the method of matched asymptotic expan- 
sions but requires a more significant numerical effort because, even though the number of 
equations to be solved at each order is the same, the operator that plays here a central role 
in the equations for the straining functions, LiipQ™^ (interpreted as an operator acting on 
f k [s,p,q}), is more complex than T>i (defined in Sect. 4.1). 



3 Surfaccs with constant s in the strained space are assumed to correspond to surfaces with constant ^ int > 
in the physical space, i.e. ip( mt "> = ip( mt \s); therefore, (with k > 0) is an ordinary differential operator 
for V> (mt) . 
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5.2. Iteration 



Thi s tech nique follows the approach taken by iPrendergast & Tomer! (Il970l ) and by 



Wilson! (119751 ). for the construction of self-consistent dynami c al mo dels of differentially ro- 
tating elliptical galaxies, and later by lLongaretti fe Lagoutd (119961 ). for their extension of 
King models to the rotating case. 



In terms of the function: 

u(r) = a[H - $ c (f)] = ^(f) + eT(f) 

inside the cluster the Poisson equation can be written as: 

9 



p(u - eT) 



while outside the cluster the Laplace equation is simply: 



V 2 u = 



(76) 



(77) 



(7f 



The boundary conditions at the origin are u(0) = ^ and Vu(0) = 0, because the tidal 
potential T(f) is a homogeneous function; the condition at large radii is u — > aH . 

The basic idea is to get an improved solution vS n ' of the Poisson equation by evaluating 
the "source term" on the right-hand side with the solution obtained in the immediately 
previous step: 



vV n > 



9 



p{u 



(n-l) 



eT) . 



(79) 



The iteration is seeded by inserting as u^°\ on the right-hand side of Eq. (I79|) . the spherical 
solution of the King models. The iteration continues until convergence is reached. 

In order to solve Eq. ( 179|) . we expand in spherical harmonics the solution and, corre- 
spondingly, the dimensionless density distribution: 



/ 



u 



in) 



r) = EE M ^mwM) 



Z=0 m=-« 
oo I 

p (n) (f) = EE£(^ 

l=Q m=-l 



Im \ 



(80) 



^1) 



so that the reduced radial problems for the functions u^{f) are: 



d 2 2d 
df 2 r df 



1(1 + 1) 



u 



(n) 
Irn 



9 -(n-l) 
Plm 



(82) 
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with boundary conditions Uqq(0) 



¥, ttJJ(O) = 0and ^(0) 



contrast with the structure of the governing equations for ip^ of Sect. 4.1 and Sect. 4.2, 
the radial part of the Laplacian appears with no "shift" , for which the homogeneous solutions 
are known analytically. Thus the full solution to Eq. (1521) can be obtained in integral form 
by the standard method of variation of the arbitrary constants: 



(«)' 

lm 



(0) = 0. Here, 



in 



u 



(*), 
00 * 



9 



r Poo 



oo '(r)dr - t / r p 00 ; (r)rf 



53) 



(2/ + l)p(tt) l 



r l I f*- l pt 1 \r')dr' + 



1 



(n ' 1) 1 



r')dr' 



14) 



The complete calculation can be found in the Appendix of iPrendergast fc Tomerl (jl97(J). 
Here we only remark that this integral form is valid in both Poisson and Laplace domains 
because it contains simultaneously the regular and the singular homogeneous solutions of 
the Laplacian. In the derivation, all the boundary conditions have been used; in particular, 
the two conditions at the origin are sufficient to obtain expression (1831) . while for expression 
the one concerning the radial derivative at the origin is used together with the one 

for I > 1). Furthermore, from the 



u 



(n) 
I in 



r(n) 



that describes the behavior at large radii (i.e. 

condition at large radii evaluated for the harmonic I = 0, i.e. w q /V4tt — > aH^ 1 ' (here the 
notation reminds us that the value of H is known only approximately and it changes slightly 
at every iteration), we find: 



Off W = 



9 



47r V47rp(vL) 



r Poo {r)dr , 



(85) 



where we should recall that beyond a certain radius p g ^ vanishes. 

In terms of the function u, the boundary of the cluster is given implicitly by: u(f) = 
eT(f). Therefore, the radial location at which the ^ vanishes is determined numerically 
from: 



2- 



pt^-^f , 9, <P) - eT(f, 9, <f>))Y lm (e, (f>)d(cose)d<P . 



16) 



-1 



In practice, to perform the iteration, the definition of a grid in spherical coordinates and 
of a suitable algorithm, in order to perform the expansion and the resummation in spherical 
harmonics of u and p, is required; the number of angular points of the grid and the maximum 
harmonic indices (/, m) admitted in the series flHUj) and (|81j) are obviously related. 
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6. Conclusions 

Spherical King models are physically justified models of quasi- relaxed stellar systems 
with a truncation radius argued to "summarize" the action of an external tidal field. Such 
simple models have had great success in representing the structure and dynamics of globular 
clusters, even though the presence of the tidal field is actually ignored. Motivated by these 
considerations and by the recent major progress in the observations of globular clusters, 
in this paper we have developed a systematic procedure to construct self-consistent non- 
spherical models of quasi-relaxed stellar systems, with special attention to models for which 
the non-spherical shape is due to the presence of external tides. 

The procedure developed in this paper starts from a distribution function identified by 
replacing, in a reference spherical model, the single star energy with the relevant Jacobi 
integral, thus guaranteeing that the collisionless Boltzmann equation is satisfied. Then 
the models are constructed by solving the Poisson equation, an elliptic partial differential 
equation with free boundary. The procedure is very general and can lead to the construction 
of several families of non-spherical equilibrium models. In particular, we have obtained the 
following results: 

• We have constructed models of quasi-relaxed triaxial stellar systems in which the shape 
is due to the presence of external tides; these models reduce to the standard spherical 
King models when the tidal field is absent. 

• For these models we have outlined the general properties of the relevant parameter 
space; in a separate paper (Varri & Bertin, in preparation) we will provide a thorough 
description of this two-parameter family of models, also in terms of projected quantities, 
as appropriate for comparisons with the observations. 

• We have given a full, explicit solution to two orders in the tidal strength parameter, 
based on the method of matched asymptotic expansions; by comparison with studies 
of analogous problems in the theory of rotating polytropic stars, this method appears 
to be most satisfactory. 

• We have also discussed two alternative methods of solution, one of which is based on 
iteration seeded by the spherical solution; together with the use of dedicated iV-body 
simulations, the ability to solve such a complex mathematical problem in different ways 
will allow us to test the quality of the solutions in great detail. 

• By suitable change of notation and physical re-interpretation, the procedure developed 
in this paper can be applied to the construction of non-spherical quasi-relaxed stellar 
systems flattened by rotation (see Appendix B). 
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• The same procedure can also be applied to extend to the triaxial case other isotropic 
truncated models (such as low-n polytropes), that is models that do not reduce to 
King models in the absence of external tides (see Appendix C). 

We hope that this contribution, in addition to extending the class of self-consistent mod- 
els of interest in stellar dynamics, will be the basis for the development of simple quantitative 
tools to investigate whether the observed shape of globular clusters is primarily determined 
by internal rotation, by external tides, or by pressure anisotropy. 



A. Details of the solution in terms of matched asymptotic expansion 

A.l. The general equation 

From the Taylor expansion about e = of the right-hand side of Eq. (TlTI) , the structure 
of the equations for ip^^ ( wr th k > 2) can be expressed as: 



V 2 + i?i(f;^) 



5> 

0=2 



r; m X, 



(Al) 



where X^j denotes the terms that arise from the derivatives of ?/>( m *)(f; e) with respect to e, 
thus expressed as products of ipi (with i = 1, .., k — 1). For fixed k and j, the quantity X^j 
is thus a sum of products of ^ m *-* with subscripts that are j-part partitions of the integer 
k. Each product of ip^""^ is multiplied by a numerical factor defined as the ratio between 
k\ and the factorials of the integers that are parts of the associated partition (if an integer 
appears m times in the partition, the factor must also be divided by ml). In particular, for 
k = 3 we have: 

X 3 , 2 = 34 nt) 4 nt) and X 3i3 = (4 mt) ) 3 , (A2) 



-3,3 



because the 2-part partition of 3 is 2 + 1 and the 3-part partition is trivially 1 + 1 + 1, thus 
the relevant equation is: 



V 2 + i?i(f;^) 



, (int) 



(int)-, 3 



(A3) 



Therefore, this formulation of the right-hand side of the general equation (tog ether with the 
term Ri(f; ^>)%p^ nt ^ on the left-hand side) brings in the Fad di Bruno formula ( Faa di Bruno 
18551 ) for the fc-th order derivative of a composite function in which the inner one is expressed 
as a series in the variable with respect to which the derivation is performed. 
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A. 2. The equation for the second order radial problem 

The expansion in spherical harmonics of (ipi^) 2 , which involves the product of two 
spherical harmonics with / = or I = 2 (with m positive and even), can be performed by 
means of the so-called 3-j Wigner symbols □. Equation fl33l) thus corresponds to the following 
set of six equations: 



T1 / ( int ) 



ti , tint) 
^2^2,20 



TI / ( int ) 

^2^2,22 



T-l / ( int ) 

^4^2,40 



ti / (int) 

v m A 2 



ti / [int) 
^4^2,44 



-R 2 (r; vl/)-— [(^«?) 2 + («) 2 + («) 2 ] 



(irai)\2 



(int) n 2n 



2V^ 



-R 2 (r',V) 



/gri,oo ri,20 lri,20 J lri,22 J 

I (int) J (int) 2V5 , 
"ri.oo "ri,22 ^~Vi,20 "ri,22 



3(^s?) 2 +i(^:22 j ) 



(int) \ 2 



70F 



(int) I (int) 

22 ) 



(A4) 
(A5) 
(A6) 
(A7) 
(A8) 
(A9) 



A. 3. The asymptotic matching for the first order solution 

To derive the first order solution, the matching between the pairs 

(^(™*),<0(M) anc i 

(^(ext) ^(iay)^ re q U j res that the internal (external) solution is expanded in a Taylor series 
around f = f tr up to terms of order 0((f — f tr ) 2 ), expressed with scaled variables, expanded 
up to 0(e 2 ), and re-expressed with non-scaled variables: 

[4 ] ] m {r,d,cj>) = 4 \r tr ) -4 Y (ftr)(ftr-f) + ^0 '"for) far ~ ff 

e + ^4\r tr ,e,4>)e 2 . (A10) 

Here the closed parentheses include either "int" or"ext", to denote internal or external 
solution, while the notation [ ]( 2 ) on the left-hand side indicates that a second order expansion 
in e has been performed. 



+ 



( ) ^ 

4 \f t r, OA) Qf-{rtr,0,<t>)(rtr ~ f) 



4 For the definition of 3-j Wigner symbols and the expression of the harmonic expansion of the product 
of two spherical harmonics, see, e.g., lEdmondsl (|l960h . Eqs. (3.7.3) and (4.6.5), respectively. 
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The boundary layer solution in the vicinity of t] = up to 0(r] 2 ), expressed with non- 
scaled variables and expanded (formally) up to third order in e is given by: 

h/^P(f , 9, 0) = F {0, 0)(f tr - f) + ^^-(hr - rf 

•tr 

+ [G (9, 0) + F 1 {8, <f>){r tr - r)]e + G x {9, 0)e 2 . 
By equating equal powers of e and (f tr — f) in flAlOj) and flAllj) . we find: 
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where the equalities on the left-hand side arise from the matching of the pair (ip( mt \ ^( /a f)) 
and identify the free angular functions (|54p - (|57p . while those on the right-hand side arise 
from the matching of the pair {ip^ xt \ip^ lav ^). We also note that (1A12I) is consistent with 
the definition of the truncation radius and that (1A14I) is equivalent to (1A13I) . because from 

Eq. m we have: rf^Vtr) = -(2/f tr ^ mt) ' (r tr ) ■ 

The free constants (l58l) -( l6T|) are thus easily determined. For a given harmonic (l,m) of 
Eqs. (1A15P and (lA16h . with / > 1, the constants A\ m and a\ m are governed by the linear 
system with i,j = 1, 2: 

M ijUj = Vi , (A18) 

where the matrix M is given by 



M 



llVtr) 



tr 



-ilihrYtr (I + I) f tr 



(l+l) 



(A19) 
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and the vectors are defined as (1*1,^2) = (Ai m , ai m ) and (t>i,t>2) = — Ti m (f tr )(l, — 2). Such 
linear system ()A18j) is well posed, i.e. detM ^ 0. 

To show this, we may integrate Eq (T3TT) for the regular solution, under the conditions 



In the vicinity of f = the quantity R\(r\ \l/)f 2 is vanishingly small, so that, if the quantity 
Ri(f;ty)f 2 remains smaller than 1(1 + 1), then the regular solution, starting positive and 
monotonic, remains a positive and monotonically increasing function of f. Indeed, we have 
checked that this condition occurs for I > 2, because for \l> 6 [0.5, 10] the quantity Ri(f; \I/)f 2 
has a maximum value in the range [4.229, 3.326]. Under this condition the function fii(f) = 
j'i(r)r + (I + l)ji(f) cannot change sign, so that detM = /^(rv) cannot vanish. This 

argument does not work for the case / = 1, in which the function fii(f) does change sign at 
a point < f < r tr , but we have checked directly that the property detM 7^ is satisfied 
also in this case. 

Then expressions ( 1521 and (1B"3"j) are easily recovered and we conclude that for the har- 
monics that are not "driven" by the tidal potential T(f) the related ai m and A\ m must vanish. 

Similarly, for a given harmonic (l,m) with I > 1, Eq. (|64p . which is obtained from 
the second order matching, and Eq. (I A 171) can be cast in the form of Eq. (IA18j) . with 
Oi,u 2 ) = (B lm ,bi m ) and («i,v 2 ) = (-gim(rtr),r tr g' lm (r tr )). Therefore, the argument provided 
above applies and we can conclude that for those harmonics for which the particular solutions 
gim(j) are absent (or, equivalently, Eq. [33] is homogeneous), the constants B\ m and bi m must 
vanish. A linear system equivalent to (1A18I) can be written for a fixed harmonic (/, m) with 
/ > 1 of the solution of general order k and the same argument applies. Therefore, we 
conclude that the k-th order term of the solution contributes only to those harmonics for 
which the particular solutions are present. 



Because we have noted (see argument introduced about the system |A18j ) that the k-th 
order term of the solution has non-vanishing contribution only in correspondence of those 
harmonics for which the component of Eq. (1A1I) is non-trivial, the discussion about the 
structure of the term reduces to the analysis of the structure of the expansion in spherical 
harmonics of the right-hand side of that equation. Recalling that the harmonic expansion of 
the product of two spherical harmonics (h,mi) and (1%, m<i) can be expressed by means of 3-j 



7l(0)=7K0) = 0: 




(A20) 



A. 4. The structure of A>th order term 
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Wigner symbols (see lEdmondslll960l . Eq. [4.6.5]), we note that the composed harmonic (l,m) 
must satisfy the following selection rules: (i) |/i — Z 2 | < I < ^i + h ("triangular inequality"), 
(ii) mi + m 2 — m i and (hi) h + h + 1 must be even. The last condition holds because in the 
cited expression the composed harmonic appears multiplied by the special case of the Wigner 
symbol with (li, I2, 1) as hrst row and (0, 0, 0) as second row. Bearing in mind that the first 
order term is characterized by harmonics with / = 0, 2 and corresponding positive and even 
values of m and that the structure of the right-hand side of Eq. (lAip can be interpreted by 
means of the partitions of the integer k, it can be proved by induction that the k-th order 
term is characterized by harmonics with I — 0, 2, 2k and corresponding positive and even 
values of m. 



B. Extension to the presence of internal rotation 

It is well known that in the presence of finite total angular momentum of the system, 
relaxation leads to solid-body rotation (e.g., see Landau & Lifchitz 1967). If we denote by uj 
the angular velocity of such rigid rotation and assume that it takes place around the z axis, 
in the statistical mechanical argument that leads to the derivation of the Maxwell-Boltzmann 
distribution one finds that in the final distribution function the single particle energy E is re- 
placed by the quantity E—ujJ z . Following this picture, we may consider the extension of King 
models to the case of internal rigid rotation. This extension is conceptually simpler than that 
addressed in the main text of this paper, because the perturbation associated with internal 
rotation, while breaking spherical symmetry, prese rves axial symmetry. We n ote that the 



models described below differ from those studied by lKormendy fc Anandl (119711). which were 



charac terized by a di f ferent , discontinuous truncation, and those by iPrendergast fc Tomer 



(119701 ) and by IWilsonl (119751 ) , which were characterized by a different truncation prescription 



and by differential rotation. 

The relevant physical model is that of a rigidly rotating isolated globular cluster char- 
acterized by angular velocity u) = ue z , with respect to a frame of reference with the origin in 
the center of mass of the cluster. We then introduce a second frame of reference, co-rotating 
with the cluster, in which the position vector is given by r = (x, y, z). In such rotating frame, 
the Lagrangian describing the motion of a star belonging to the cluster is given by: 

£ = i(i 2 + y 2 + z 2 + 2uyx - 2uxy) - $ cen (x, y) - $ c (x, y, z) , (Bl) 

where $ cen (x, y) = —(x 2 + y 2 )u 2 /2 is the centrifugal potential; the energy integral of the 
motion (called the Jacobi integral) is: 

H = l -{x 2 + y 2 + z 2 ) + $ cen + $ c . (B2) 
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As in the tidal case, the extension of King models is performed by considering the distribution 
function f K (H), as in (jTJ), for H < H and fx^H) = otherwise, with H the cut-off constant. 
The dimensionless energy is defined by: 

V>(r) = a{H - [$ c (r) + $ cen (x, y)]} , (B3) 

and the boundary of the cluster, implicitly defined as if)(r) — 0, is an equipotential surface 
for the total potential + $ cen - Its geometry, reflecting the properties of the centrifugal 
potential, is characterized by symmetry with respect to the z-axis and reflection symmetry 
with respect to the equatorial plane (x, y) . The constant ip family of surfaces, much like the 
Hill surfaces of the tidal case, is characterized by a critical surface which distinguishes the 
closed from the opened ones and in which the points on the equatorial plane are all saddle 
points. In these points the centrifugal force balances the self-gravity of the cluster; their 
distance from the origin, which we call break-off radius (vb), can be determined from the 
condition: 

^(r B ,0,0) = ^(0,r B ,0) = 0. (B4) 
ox ay 

Following the argument that we gave in Sect. 2.2 for the tidal radius, we find a zero-th order 
approximation for the break-off radius: = (GM/oo 2 ) 1 ^ 3 , where M is the total mass of 
the cluster. The discussion about the parameter space that characterizes these models is 
equivalent to the one presented in Sect. 2.3, with: 

(B5) 

a dimensionless parameter that measures the strength of the rotation with respect to the 
central density of the cluster, playing the role of e. Similarly, we may define an extension 
parameter a = r tr /rB (instead of S, see Sect. 2.3). For every value of the dimensionless central 
potential well a maximum value for the strength parameter, Xcr, exists, corresponding to 
the critical condition in which the boundary of the cluster is given by the critical equipotential 
surface. 

The relevant equations in the construction of these (fully) self-consistent models can be 
expressed in dimensionless form by means of the same rescaling of variables performed in 
Sect. 3. Thus, for ip > 0, the Poisson equation can be written as: 



V 2 ^ = -9 



m_ 2x 



(B6) 



where p is defined as in (11 II) . For negative values of ip we should refer to: 

V 2 ^ = 18 X , (B7) 
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with the boundary conditions at the origin written as in f|T9l) - fl20l) and at large radii given 
by ip + xC — > aH , where C = — (9/2)(x 2 + y 2 ) is the dimensionless centrifugal potential. 

The solution up to second order in terms of matched asymptotic expansions presented 
in Sect. 4 can be adapted to this case without effort. In fact, with respect to the calcula- 
tion presented in the main text only two differences occur: (i) wherever the constant term 
— V 2 T = —9(1 — v) appears, it must be replaced here by — V 2 C = 18 (the sign is the same 
in the two cases, because 1 — v < 0, see Sect. 2.1), and (ii) thanks to axisymmetry, in 
the angular part of the Laplacian the derivative with respect to the toroidal angle can be 
dropped and thus the terms of the asymptotic series (1221 - fl23|) can be expanded by means of 
Legendre polynomial^ instead of spherical harmonics. The latter property implies that the 
radial part of each term of the asymptotic series is characterized by only one index, I, i.e. 
we can write We note that the differential operator that appears on the left-hand side 
of the relevant equations for the solution defined in the internal region is still T>\, and thus 
also the functions 7/(r) can be introduced in the same way as before. As to the equations 
corresponding to (|5T]) - (!53|) . the formal solutions of the equations in the boundary layer, the 
angular functions Fj and Gi (with i — 0, 1, 2) now depend only on the poloidal angle 6. 
About the external solution, an expression analogous to (13T)j) can be used, with the particu- 
lar solution given by \C instead of eT. The centrifugal potential contributes, as in the case 
of the tidal potential, only with monopole and quadrupole terms, explicitly: 



Finally, as a result of the matching of the pair (ip ( - mt \ ^A'°tf)) up to second order, the expres- 
sions for the angular functions ( !54l) - (!57j) and (|64l) -(j65ll are still applicable. In addition, from 
the matching of the pair (ip( lay \ ip ext ) up to second order, we find that the explicit expressions 
for the free constants follow Eqs. ( I58l) - (l63l) and ( l66l) -( J7TT) . provided that we drop everywhere 
the index m and we replace 3T 00 (ft r ) / (2\/tt) with 3C (f tr ) j \[2 in fl60|) and T 00 (r tr ) / \-Jiif tr ) 
with 2Co(r tr ) J \\/2f tr ) in fl6T|) . Obviously, the particular solutions f and gi (with I = 0,2, 4) 
are different from the ones obtained in the tidal case, because the right-hand side of the 
relevant equations is different. Also in this case, it can be proved by induction that the k-th 
order term has non- vanishing contributions only for I = 0,2, 2k. 

For completeness, we record the explicit expression of the second order equations in the 



(B8) 




(B9) 



5 Following I Abramowitz fc Stegunl (|19641 ). we use Legendre polynomials as denned in (22.3.8), i.e. with 
Condon-Shortley phase, and normalized with respect to the relation (22.2.10). We remark that, although 
they are structurally equivalent to zonal spherical harmonics, the normalization is different. 
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We remark that the Legendre expansion of the product of two Legendre polynomials is 
straightforward, because the 3-j Wigner symbols of interest all belong to the special case 
with (0, 0, 0) as second row. 



C. Extension of other isotropic truncated models 

The procedure developed in Sects. 3 and 4 can be applied also to extend other isotropic 
truncated models, different from the King models, to the case of tidal distortions. Here 
we briefly describe the case of low-ra polytropes (1 < n < 5), which are particularly well 
suited to the purpose, because they are characterized by a very simple analytical expression 
for the density as a function of the potential; this class of models was also conside red by 



Weinberg] (119931 ). In the distribution function that defines the polytropes (e.g., see iBertin 
2000l ). we may thus replace the single star energy with the Jacobi integral (see definition [5]) 
and consider: 

f P (H) = A(H - H) n ^ 2 , (CI) 

for H < H , and a vanishing distribution otherwise. Unlike the King family discussed in the 
main text, these models have no dimensionless parameter to measure the concentration of 
the stellar system, which depends only on the polytropic index n; in fact, the spherical fully 
self-consistent polytropes are characterized only by two physical scales, which are associated 
with the cut-off constant H and the normalization factor A. Below we consider values 
of n < 5, so that the models have finite radius. Therefore, the relevant parameter space 
for the tidally distorted models is represented just by the tidal strength parameter e (see 
definition [13]), which, for a given value of the index n, has a (maximal) critical value. The 
definition (TBI) for the extension parameter 5 is still valid if by r tr we denote the radius of 
the unperturbed spherical configuration. The associated density functional is given by: 

pW = p r , (C2) 

where the dimensionless escape energy is given by: 

fj ' (C3) 
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with c n = (27r) 3 / 2 r(n — 1/2) A/n\. The boundary of the perturbed configuration is defined 
by ip(r) = 0, following the same arguments described in the main text. Here po can be 
interpreted as the central density if we set ip(0) = 1. 

For ip > 0, the relevant equation for the construction of the self-consistent tidally 
distorted models is the Poisson equation, which, in dimensionless form, is given by: 

v 2 ^ = -\r-<i-v)] , (C4) 

while for negative values of ip we must refer to Eq. (Tl8l) . Here the rescaling of variables has 

been performed by means of the scale length ( = J p l J n ~ l / (AnGcrl™). The relevant boundary 
conditions are given by ip(0) = 1 instead of (fT9l) . while (1201) and (T2"Tj) hold unchanged. 

If the polytropic index is in the range 1 < n < 5, the solution up to second order 
presented in Sect. 4 is fully applicable, provided that we note that the problem for the 
zero-th order term o f the series (1221) is now given by the Lane-Emden equation (see, e.g. 



Chandrasekharl 1 1 9391 ) : 



4 nt) + -A ni) = - W int) ) , (C5) 



r 

with ipQ nt \o) = 1 and ip^™^ (0) = 0, where the symbol ' denotes derivative with respect to 
the argument f; explicitly, the truncation radius f tr is now defined by ipQ nt \r tr ) = 0, i.e. it 
represents the radius of the so-called Emden sphere. Correspondingly, the quantities called 
Rj in the main text must be re-defined as: 

the value of j at which the quantity Rj may start to diverge depends on the index n. 
Obviously, in Eq. (143]) . i.e. in the Poisson equation defined in the boundary layer, p(er) 
must be replaced by {er) n . This makes it clear that the value of the polytropic index n 
directly affects the order, with respect to the perturbation parameter, at which the density 
contribution on the right-hand side of Eq. (145!) comes into play and therefore changes the 
matching procedure. If n > 1, the density contribution emerges only after the second order 
and thus the full procedure described in Sect. 4 is valid. In contrast, if n < 1 the procedure 
described in the main text is applicable only up to first order while the calculation of second 
order terms would require a re-definition of the boundary layer (as it happens for the case 
discussed in the main text when terms of order k > 3 are desired). In closing, we note 
that the procedure presented in this paper can be applied also to isotropic truncated models 
with more comp licated expressi ons for the density functional (e.g. the family of models 



fi n proposed by iDavoustl (119771 ). without boundary conditions on tangential velocity, for 
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which the density functional are expressed in terms of the error function and of the Dawson 
integral), bearing in mind the last caveat about the possibility that the density contribution 
may affect at some order the boundary layer, thus requiring a reformulation of the results 
presented in Sect. 4. 
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Fig. 1. — The critical Hill surface (in dimensionless variables) for a second-order model 
with \I> = 2 and e = 7.043 x 10~ 4 (corresponding to 5$ = 0.671); the galactic potential is 
Keplerian {v — 3). 
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Fig. 2. — Sections in the three coordinate planes for a second-order model with ^ = 2 
and e = 7.000 x 10~ 4 , characterized by high tidal distortion (5 = 0.669 ~ see Fig. [T]) 
illustrating the boundary surface of the triaxial model (solid), of the internal region (dotted), 
and of the corresponding spherical King model (dashed); the filled area represents the inner 
region of the boundary layer. Note the compression and the elongation with respect to the 
unperturbed configuration in the z- and ^-direction, respectively. The galactic potential is 
Keplerian [y — 3). 



